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Abstract 



The closest tensors of higher symmetry classes are derived in explicit form for 
a given elasticity tensor of arbitrary symmetry. The mathematical problem is to 
minimize the elastic length or distance between the given tensor and the closest 
elasticity tensor of the specified symmetry. Solutions are presented for three distance 
functions, with particular attention to the Riemannian and log-Euclidean distances. 
These yield solutions that are invariant under inversion, i.e., the same whether 
elastic stiffness or compliance are considered. The Frobenius distance function, 
which corresponds to common notions of Euclidean length, is not invariant although 
it is simple to apply using projection operators. A complete description of the 
Euclidean projection method is presented. The three metrics are considered at a 
level of detail far greater than heretofore, as we develop the general framework 
to best fit a given set of moduli onto higher elastic symmetries. The procedures 
for finding the closest elasticity tensor are illustrated by application to a set of 21 
moduli with no underlying symmetry. 

1 Introduction 

We address the question of finding the elastic moduli with a given material symmetry 
closest to an arbitrary set of elastic constants. There are several reasons for reducing a 
set of elastic constants in this way. One might desire to fit a data set to the a priori known 
symmetry of the material. Alternatively, a seismic simulation might be best understood in 
terms of a model of the earth as a layered transversely isotropic medium, even though core 
samples indicate local anisotropy of a lower symmetry. More commonly, one might simply 
want to reduce the model complexity by decreasing the number of elastic parameters. In 
each case a distance function measuring the difference between sets of elastic moduli is 
necessary to define an appropriate closest set. The most natural metric is the Euclidean 

* Laboratory for Mathematical and Numerical Modeling in Engineering Science, National Engineering 
School at Tunis, ENIT-LAMSIN, B.P. 37, 1002 Tunis Belvedere, Tunisia, maher. moakher@enit.rnu.tn 

tRut gers University, Department of Mechanical and Aerospace Engineering, 98 Brett Road, Piscat- 
away, NJ 08854-8058, norris@rutgers.edu 



MOAKHER AND NORRIS 



February 4, 2008 



2 



norm in which the "length" || C|| of a set of moduh-"^ Cijki is || C|| = (CijkiCijkiY^"^ . The Eu- 
chdean distance function is, however, not invariant under inversion of the elasticity tensor 
- one obtains a different result using compliance as compared with stiffness. Alternative 
distance functions have been proposed which do not have this failing. In particular, the 
Riemannian distance function of Moakher P and the log-Euclidean metric of Arsigny et 
al. [2j are invariant under inversion, and as such are preferred general elastic norms. 

The problem of simplifying elastic moduli by increasing the elastic symmetry was 
apparently first considered by Gazis et al. [3]. They provided a general framework for 
defining and deriving the Euclidean projection, i.e., the closest elasticity using the Eu- 
clidean norm. Several examples were given, including the closest isotropic material, which 
agrees with the isotropic approximant derived by Fedorov Fedorov obtained isotropic 
moduli using a different criterion: the isotropic material that best approximates the elas- 
tic wave velocities of the given moduli by minimizing the difference in the orientation 
averaged acoustical tensors. It may be shown j3] that the generalization of Fedorov's 
criterion to other symmetries is satisfied by the Euclidean projection of the stiffness ten- 
sor onto the elastic symmetry considered. The Euclidean projection is also equivalent to 
operating on the given elasticity tensor with the elements of the transformation group of 
the symmetry in question [Hj. This approach has been used by Frangois et al. [2j to find 
the closest moduli of trigonal and other symmetries for a set of ultrasonically measured 
stiffnesses. The equivalence of the projection and group transformation averaging will 
be discussed in Section El The Euclidean projection approach has received attention in 
the geophysical community from modelers interested in fitting rock data to particular 
elastic symmetries |HliniE3^3- The work of Helbig jH] is particularly comprehensive. 
He considers the problem both in terms of the 6x6 matrix notation [T^ IT!^ and in terms 
of 21- dimensional vectors representing the elastic moduli. The latter approach has been 
developed further by Browaeys and Chevrot jT3] who describe the projection operators 
for different elastic symmetries in the 21-dimensional viewpoint. Bellinger ^H] presents 
an algorithm for finding the closest transversely isotropic medium by searching over all 
orientations of the symmetry axis. Bellinger et al. ^H] proposed a distinct approach to 
the problem of finding closest elasticity, based on the idea of generalized rotation of the 
21-dimensional elasticity vector. 

Although the Euclidean projection is commonly used in applications, it suffers from 
the fundamental drawback alluded to earlier, i.e., it is not invariant under inversion of 
the elasticity tensor. The projection found using the elastic stiffness is different from that 
obtained using the compliance tensor. This fundamental inconsistency arises from the 
dual physical properties of the elasticity tensor/matrix and its inverse. The projection of 
one is clearly not the same as the inverse of the projection of the other, although both 
projections by definition possess the same elastic symmetry. While there are circumstances 
in which the Euclidean projection is preferred, e.g., for the generalized Fedorov problem 
of finding the best acoustical approximant of a given symmetry jH], there is a clear need 
for a consistent technique to define a "closest" elastic material of a given symmetry. 
The solution to this quandary is to use an elastic distance function that is invariant 
under inversion. Several have been proposed of which we focus on two, the Riemannian 
distance function due to Moakher and the log- Euclidean length of Arsigny et al. [2], 
described in Sectional Moakher P describes how the Riemannian distance function gives 
a consistent method for averaging elasticity tensors. The only other application so far of 
these invariant distance functions to elasticity is by Norris who discusses the closest 



""^ Lower case Latin suffices take on tfie values 1, 2, and 3, and the summation convention on repeated 
indices is assumed. 
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Figure 1: The sequence of increasing elastic symmetries, from triclinic to isotropic. 

The number of independent elastic constants are listed. The dashed boxes 
for trigonal 6 and tetragonal 6 indicate that these are obtained from the 
lower symmetries by rotation, and do not represent new elastic symmetries. 
The true symmetries are the eight in solid boxes |18j . 



isotropic elasticity. 

The outline of the paper is as follows. The three distance functions are introduced in 
Section El after a brief review of notation in Section El Section jH defines the different types 
of elastic symmetry using the algebraic tensor decomposition of Walpole ^Hj- Euclidean 
projection is presented in Sectional with results for particular elastic symmetries summa- 
rized in Appendix El Before discussing the closest tensors using the logarithmic norms, 
methods to evaluate the exponential and other functions of elasticity tensors are first de- 
scribed in Section ini Section (3 discusses the closest tensors using the logarithmic norms, 
with particular attention given to isotropy and cubic symmetry as the target symmetries. 
Numerical examples are given in Section |H1 and final conclusions in Section (HI 



2 Notation 

Elasticity tensors relate stress T and infinitesimal strain E linearly according to 

T = CE, E = §T, (1) 

where C and S denote the fourth-order stiffness and compliance tensors, respectively. 
They satisfy C§ = §C = I, the identity. Although we are concerned primarily with 
fourth-order tensors in 3-dimensional space, calculation and presentation are sometimes 
better performed using second-order tensors in 6-dimensions. Accordingly we define 
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S{d, r) as the space of symmetric tensors of order r in (i— dimensions. Elasticity ten- 
sors, denoted by Ela C iS(3,4), are positive definite, i.e., A G Ela if (B, AB) > for 
all nonzero B G 5(3,2). Components are defined relative to the basis triad {01,62,63}; 
thus, a = ajGj, A = Aij 6j ® e^-, and A = Aij^i 6j ® e^- ® 6^ ® e^, where the summation 
convention is assumed over 1,2,3 for lower case subscripts. Symmetry of second-order 
(r = 2) tensors implies Aij = Aji, while for A G 5(3, 4) the elements satisfy 

The first pair of identities reflects the symmetry of the stress and the strain, while the last 
one is a consequence of the assumed existence of a strain energy function, and consequently 
elasticity tensors have at most 21 independent components. 

Throughout the paper lower case Latin, upper case Latin and "ghostscript" indicate 
respectively 3— dimensional vectors and tensors of order 2 and 4; e.g., vector b, A G 
iS(3,2), A G iS(3,4). The basis vectors are assumed orthonormal, 6j ■ ej = 6ij, so that 
products of tensors are defined by summation over pairs of indices: (AB)jj = AikBkj and 
{KB>)ijki = AijpqBpqki- The inner product for tensors is defined as 

{u, v) = tT{uv), (3) 

where tr A = An, and for elasticity tensors, tr A = Aijij. The norm of a tensor is 

\\u\\ ^ {u, uy/\ (4) 



We take advantage of the well known isomorphisms between 5(3,2) and 5(6, 1), and 
between 5(3,4) and 5(6,2). Thus, fourth-order elasticity tensors in 3 dimensions are 
equivalent to second-order symmetric tensor of 6 dimensions [20] , with properties B A ^ 
Ba, AB ^ AB, and ( A, B) = (A, B). Vectors and second-order tensors in six dimension 
are distinguished by a hat, e.g., vector a, A G 5(6, 2). Components are defined relative to 
the orthonormal sextet {sj, / = 1, 2, . . . , 6}, 6/ ■ ej = by a = ajSj, A = Ajjej ® ej, 

with the summation convention over 1,2, ... ,6 for capital subscripts. Also, tr A = Ajj. 
The connection between 5(6, 1) and 5(3, 2) is made concrete by relating the basis vectors: 

6 - vector dyadic 

ei 64 ^ ei(g)ei ^(62(8)63 + 63062) (5) 

62 65 62 62 ^(63(8)61+61(8)63) 

63 66 63(863 ^(61(862 + 62(861) 

This implies a unique A G 5(6,2) for each A G 5(3,4), and vice versa. Let C G Ela be 
the tensor of elastic stiffness, usually defined by the Voigt notation: Cijki = cjj, where / 
or J = 1, 2, 3, 4, 5, 6 correspond to ij or kl = 11, 22, 33, 23, 13, 12, respectively, and c/j are 
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the elastic moduli in the Voigt notation 
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(6) 



The isomorphism implies that the associated matrix C G 5(6,2) is positive definite and 
has elements 



/ Cll 


C12 


Cl3 


22C14 


2^Ci5 


25Ci6\ 
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(7) 



The spectral decomposition of elasticity tensors can be expressed in both the fourth- 
order and second-order notations [1211201, 



1=1 



1=1 



where A/ and hj, I = 1,2, ... ,6 are the eigenvalues and eigenvectors of the matrix C, N/ 
are the associated dyadics, with {hj, hj) = (Nj, Nj) = 6ij. Also, A/ > by virtue of the 
positive definite nature of the strain energy. 

The elastic moduli C can also be expressed as a vector with 21 elements 



X = {xi,X2,X3, X4, X5, X6, X7, X8, Xg, Xio , 2;ii , X12 , X13, X14, Xi5,XiQ,Xn, Xis, X19, j;20, X21Y 
= (cii,C22, C33 , \/2c23 , \/2ci3 , V2ci2 , 2C44 , 2C55 , 2cm , 2ci4 , 2C25 , 2C36 , 

2C34, 2ci5, 2C26, 2C24, 2C35, 2ci6, 2\/2c56, 2\/2c46, 2\/2c45^ 

= (cii , C22 , C33 , \/2c23 , \/2ci3 , \/2ci2 , C44 , C55 , C66 , V^Cu , \/2c25 , \/2c36 , 

\/2c34 , \/2ci5 , V2C26 , V2C24 , V2cs5 , V2ciQ , \/2c56 , V2Ci6 , \/2c45') • 



(9) 
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The V2's ensure that the inner product preserves the norm of the elastic moduh, whether 
it is a tensor, matrix or vector. Thus, the norm of an elasticity tensor can be expressed 
in a various ways depending on how C is represented, 

II C||' = ( C, C) = Q,kiC,,M = cjjcij = AjAj = \\Xf . (10) 

3 Elastic distance functions 

We consider three metrics for Ela: the Euclidean or Frobenius metric dp, the log- 
Euclidean norm d^ j2] and the Riemannian metric da pQ. They are defined for any 
pair of elasticity tensors as 

d^(Ci, C2) = ||Ci - C2II, (lla) 
dL(Ci, C2) = ||Log(Ci) -Log(C2)|| , (lib) 



dK(Ci, C2) 



Log(C^^/^ C2 C7^/^) 



'llcl 



The Riemannian distance d^ is related to the exponential map induced by the scalar 
product on the tangent space to Ela at C [TJ The log- Euclidean metric dL can be 
motivated by the following definition of tensor multiplication |2] : 

Ci0 C2 = exp(Log(Ci) +Log(C2)) . (12) 

The product preserves symmetry and positive definiteness, unlike normal multiplica- 
tion. The three metrics have the usual attributes of a distance function d: 

1. non-negative, d(Ci, C2) > with equality iff Ci = C2, 

2. symmetric in the arguments, d(Ci, C2) = d(C2, Ci), 

3. invariant under a change of basis, d(C|, Cg) = d(Ci, C2) for all proper orthogonal 
coordinate transformations {61,62,63} {61,62,63}, and 

4. it satisfies the triangle inequality, d(Ci, C3) < d(Ci, C2) + d(C2, C3). 

The Riemannian and log-Euclidean distances possess the additional properties that 
they are invariant under inversion and (positive) scalar multiplication 

dL,ii(Cr\ C2') =di,^(Ci, C2), and di,^(a C, a C) = dz.,^( Ci, C2), a > 0. (13) 

The following inequality between the logarithmic distance functions is a consequence of 
the metric increasing property of the exponential |2Hj 



dL(Ci, C2) <dK(Ci, C2). (14) 
We also have (see [23] for the inequality) 

di(C?, C^) = |6|dL(Ci, C2), beR (15a) 

dBi Cl C^) < b d«( Ci, C2) , be [0, 1], (15b) 
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and the Riemannian distance is invariant under congruent transformations, i.e., 

dR(TCiT*, TC2T*) = dij(Ci, C2), VT invertible. (16) 

Distance functions satisfying (fTSj) are called bi-invariant, a property that makes the 
Riemannian and log-Euchdean distances consistent metrics for elasticity tensors. The 
Frobenius norm, not being invariant under inversion, gives different results depending 
on whether the stiffness or compliance is considered. Other metric functions may be 
considered. For instance, the KuUback-Leibler metric PP is invariant under inversion and 
congruent transformation but does not satisfy the triangle inequality, which we take here 
as prerequisite for consideration as a distance function. 

Each distance function has a geometrical interpretation. For instance, the midpoint 
between Ci and C2 is defined as the unique C3 such that d(Ci, C3) = d(C3, C2) = 
id( Ci, C2). Using dp, the midpoint is the vector halfway between the 21— vectors Xi and 

X2 as defined by ©, i.e., C3 = i( Ci+ C2). The midpoint using d^ is C3 = ( Ci0 €2)^/^. 
The midpoint with d^j is j21] C3 = Ci( C^^^ C2)^^^ which can be expressed in several other 
ways, see [T]. More generally, the midpoint is the value of the geodesic C{t) at t = ^, 
where 

{(1 - t) Ci + t C2, Frobenius, 
exp((l -t)Log(Ci) +tLog(C2)), log - Euclidean, < t < 1. (17) 
Ci exp (t Log( Cj^^ C2)) , Riemannian, 

4 Elasticity tensors of the different symmetry classes 

Elasticity tensors for the different symmetry classes are described in this section. The 
choice of representation of Ela is important for discussing the closest approximants, par- 
ticularly for the bi-invariant metrics. It should preferably be independent of the coordinate 
system, although at the same time, practical application is normally in terms of the Voigt 
notation, so the format should not deviate too far from this. We begin with a review of 
alternative representations for elastic moduli. 

4.1 Representation of elasticity tensors 

There are many representations of Ela in addition to the standard Voigt notation of eq. 
(jH)). The 21— dimensional vector format ^Ij is useful for some applications, including the 
Euchdean projection, as we show later. However, these two representations depend upon 
the coordinate system. Among the coordinate-free forms of Ela that can be identified in 
the literature we distinguish (i) spectral decomposition, (ii) algebraic decomposition, (iii) 
groups and reflection symmetries, (iv) harmonic decomposition, (v) integrity bases. 

The first, (i), spectral decomposition dates back to Kelvin [2^]. The idea is simple: 
the elasticity tensor operates on the six-dimensional space of symmetric second-order 
tensors, and therefore has a six-dimensional spectral form. The associated eigenvalues 
are called the Kelvin moduli [T21II31- Recent related developments are due to Rychlewski 
[T^ and Mehrabadi and Cowin [201 El, who independently rediscovered and extended 
Kelvin's approach, see also ||26ji2Z]. The associated six-dimensional tensor representation 
will be used here for practical implementation, (ii) Tensor functions operating on Ela 
can be greatly simplified using the irreducible tensor algebra proposed by Walpole jTHI, 
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and independently by Kunin j2H]- This procedure provides the most efficient means of 
representing elastic tensors of a given symmetry class, especially those of high symmetry. 
We adopt this method to develop most of the results here, (iii) The group of rotations 
associated with elastic symmetry provides an irreducible representation [21]. There are 
various related ways of considering elasticity tensors in terms of rotational group properties 
of tensors, e.g., based on Cartan decomposition jHOI, complex vectors and tensors jSHEZI, 
and subgroups of 0(3) [SHI El! • These ideas are closely related to definitions of elastic 
symmetry in terms of a single symmetry element: refiection about a plane. The necessary 
algebra and the relationship to the more conventional crystallographic symmetry elements 
are described in detail by Cowin and Mehrabadi [SH]- Comparison between the rotation- 
based approach and that using symmetry planes is provided by Chadwick et al. jl8j . 
(iv) Backus proposed a representation of Ela in terms of harmonic tensors. These 
are based on an isomorphism between the space of homogeneous harmonic polynomials of 
degree q and the space of totally symmetric tensors of order q. There has been considerable 
interest in Backus's representation [SZl EHl lEl IHO]- For instance, Baerheim and Helbig 
|in] provide an orthonormal decomposition of Ela in terms of harmonic tensors, (v) 
Elasticity decomposition via integrity bases has been studied considerably [HI HJl 
El 1^- An integrity basis is a set of polynomials, each invariant under the group of 
symmetry transformations, such that any polynomial function invariant under the group 
is expressible as a polynomial in elements of the integrity basis. For instance, Tu used 
an integrity basis to construct five hierarchies of orthonormal tensor bases which span the 
space of elastic constants of all crystal systems. Any elastic tensor of order four possessing 
certain symmetry may be decomposed into a sum of tensors of increasing symmetry. The 
resulting decomposition has considerable similarity to the decomposition generated in 
Section El using Walpole's irreducible elements with Euclidean projection. Among other 
methods for decomposing elasticity tensors, we mention the scheme of Elata and Rubin 
|in] who use a set of six vectors related to a regular icosahedron to form a basis for Ela. 
This basis naturally splits a tensor into deviatoric and non-deviatoric parts. 

We use Walpole's fH] algebraic representation as it provides a consistent and straight- 
forward means to define projections of fourth-order tensors onto the given elastic sym- 
metry. This representation is better for our purposes than the spectral decomposition, 
since it is independent of the elastic moduli (i.e., no distributors ^H]) and depends only 
on the crystallographic orientation. It uses the notion of basis tensors, similar to but not 
the same as for a vector space, which makes it very suitable for Euclidean projection. We 
work mainly with fourth-order tensors directly, although the associated six-dimensional 
matrix notation is also provided. The latter is simpler for purposes of computation, e.g., 
many of the matrix operations are easily implemented using MATLAB. We begin with 
the highest elastic symmetry. 

4.2 Isotropic system 

A general isotropic fourth-order tensor is given by 



A = aJ + 6K, 



(18) 



where J and K are the two linearly independent symmetric tensors defined by 



J = 



K= I- J. 



(19) 
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The component forms follow from lijki = ^{SikSji + SuSjk) and Jijki = ^SijSki- Note that 
the tensors J and K sum to the identity, 

I = J + K, (20) 

and they satisfy the multiplication table 

f = J, = K, JK = KJ = O. (21) 

The Euclidean lengths are || J|| =1, || K|| = -\/5. Tensors J and K form, respectively, one 
and five- dimensional Kelvin subspaces [13] . 

An important result (see [SI E]) is that the closest isotropic elasticity tensor in the 
Euclidean sense is 

Ciso = 3kJ + 2^K, (22) 

where 

9k = 3( C, J) = Cii + C22 + C33 + 2Ci2 + 2Ci3 + 2C23, ^^^^ 

30/i = 3( C, K) = 2(cii + C22 + C33 - C12 - C23 - C31) + 3(c44 + C55 + cee) • 

Gazis et al. jS| obtained ()22|) using methods similar to those we will generalize in Sectional 
Fedorov actually found k and /i using an apparently different approach - the minimization 
of the mean square difference in the acoustical or Christoffel matrices for the original and 
isotropic systems. However, it can be shown j3] that the generalization of Fedorov's 
approach to other symmetries is identical to the Euclidean minimization. The moduli 
will be derived later in the context of the general theory for Euclidean projection. 



4.3 Cubic system 

Let a, b and c be three mutually orthogonal unit vectors that describe the three crystal- 
lographic directions of a cubic medium. We introduce the second-order tensors 

U = ^(a® c + c®a), V = ^(bOc + cOb), W = ^(a O b + b ® a), 

X = ^(cOc-a® a), Y = ^(b ® b - c ® c), Z = ^(a (g) a - b O b), ^'^^^ 

and fourth-order tensors L and M defined by them, 

L = U®U + V(8)V + W®W, M = -(X(8)X + Y(g) Y + Z® Z). (25) 

3 

The tensors J, L and M sum up to the identity tensor, and L and M partition K, i.e., 

I=J+L+M, K=L+M. (26) 
The multiplication table for cubic tensors is 

f = J, = L, = M, JL = LJ = JM = MJ = LM = ML = O, (27) 

and the Euclidean lengths are ||L|| = V3, ||M|| = ^2. 

A general symmetric fourth-order tensor for a cubic media is given by a linear combi- 
nation of the three linearly independent tensors J, L and M 

A = aJ + 6L + cM. (28) 

The tensors J, L and M form, respectively, one, three and two-dimensional Kelvin sub- 
spaces [13]. 
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4.4 Transversely isotropic system 

Assume that the unit vector c characterizes the preferred direction of a transversely 
isotropic medium. Let P and Q be the second-order tensors 



1 



P = c(g)c, Q = — (I-c®c), 
and define the six hnearly independent fourth-order elementary tensors 



(29) 



El 
F = 



= P 
W 



w + z$ 



Q ® Q, E3 



= P ® Q, E4 = Q 
U O U + V O V. 



(30) 



Both El and E2 are symmetric but the tensors E3 and E4 are not, although they are 
the transpose of one another and so their sum is a symmetric tensor. These asymmetric 
tensors are introduced because the set {Ei, E2, E3 -|- E4} is not closed under tensor 
multiplication. The multiplication table for the six elementary tensors is 





El 


E2 


E3 


E4 


F G 


El 


El 
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E2 





E2 





E4 





E3 





E3 





El 





E4 


E4 





E2 








F 
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G 














G 



(31) 



The tensors Ei, E2, F and G sum up to the identity, 

I = El + E2 + F + C 



(32) 



The Euclidean lengths are || Ei|| = || E2II = 1, || E3 + E4II = || F|| = || G|| = V2. 

A general symmetric fourth-order tensor of transversely isotropic symmetry is given 
by the 5-parameter linear combination 



aEi + 6E2 + c(E3 



E4 



f¥ + gG. 



(33) 



Note that the base tensors do not correspond to the Kelvin modes, although some can 
be identified as such, e.g., G is a two-dimensional Kelvin subspace The three bases 
tensors {Ei, E2, E3 -|- E4} together define a two-dimensional subspace. 



4.5 Tetragonal system 

Nine base tensors are required: {Ei, E2, E3, E4, Fi, F2, F3, F4, G}, where the new 
tensors are a symmetric pair Fi, F2, and a pair that are mutual transposes, F3, F4, 

Fi = W®W, F2 = Z®Z, F3 = W®Z, F4 = Z®W. (34) 

In comparison with the transversely isotropic system, F is replaced by the symmetric pair 
of tensors, i.e., F = Fi -|- F2, and the decomposition of the identity is 



I = El + E2 + Fi + F2 + G. 



(35) 
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The multiplication table of the nine elementary tensors is 





El 


E2 


E3 


E4 


± 


F2 


F3 


F4 


G 


El 


El 
























E2 





E2 




E4 

















E3 





E3 




El 

















E4 


E4 





E2 




















Fi 













Fi 





F3 








F2 
















F2 





F4 





F3 
















F3 





Fi 





F4 













F4 





F2 








G 


























G 



(36) 



Under multiplication the Ej's and Fj's decouple from one another and from G. Further- 
more, the algebra of the Fj's is similar to that of the Ej's. The Euclidean lengths of the 
new tensors are || Fi|| = || F2II = 1, || F3 + F4II = V^. 

A general symmetric fourth-order tensor of tetragonal symmetry is given by the 
6— parameter combination 

A = aEi + 6E2 + c(E3+ E4) +pFi + gF2 + r(F3 + F4)+^G. (37) 

Tetragonal symmetry is often represented by 6 rather than 7 independent parameters. 
Fedorov j3] pointed out that transformation by rotation about the c axis by angle 6 yields 
zero for the transformed coefficient r if (see eq. (9.7) in |3]) 

tanAO = r/{q-p). (38) 

The transformation (jHH|) depends on knowledge of the coefficients. However, it is assumed 
here that the only properties of the tetragonal symmetry available a priori are the orthog- 
onal crystal axes a, b and c. The reduction from 7 to 6 parameters can be achieved after 
the effective tetragonal material is found, but the axes of the final 6-parameter material 
depend upon the initial moduli. Hence, it is important to retain the seven-parameter 
representation (jHIj). 



4.6 Trigonal system 

Trigonal symmetry is characterized by three planes of reflection symmetry with normals 
coplanar and at 120° to one another. This symmetry class is also known as hexagonal fH] 
but we refer to it as trigonal since hexagonal symmetry is often used synonymously with 
transverse isotropy. It is advantageous to choose non-orthogonal basis vectors a', b', c 
with a' and b' normals to two of the planes and c perpendicular to them (the pair a' 
and b' may also be represented in terms of the orthonormal vectors via, for instance, 
a' = a, b' = — |a -|- ^b). In order to describe symmetric fourth-order tensors of trigonal 
symmetry, following ^H] we introduce the second-order tensors 



'|(a'®b' + b'®a' + a'®a'), T = y'Ka' ® b' + b' ® a' + b' ® b'), (39a) 
U' = ^(a'® c + c®a'), V = -^(b'®c + c®b'), (39b) 
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the two symmetric fourth-order tensors 

Ml = |(S®S + T®T- iS®T- iT® S), (40a) 
M2 = |(U' ® U' + V ® V - |U' ® V - |V' ® U') , (40b) 
and two pairs of mutually transpose tensors: 

M3 = |(S ® U' + T ® V - It ® U' - |S ® V) , (41a) 
M4 = |(U' ® S + V ® T - iU' ® T - iV ® S) , (41b) 

and 

M5 = -|(S® V'-T®U'), M6 = ;|(V'®S-U'®T). (41c) 

To make the Mj's closed under multiplication, we need to further introduce two skew- 
symmetric tensors 

M7 = -|(T®S-S®T), Rg = ^(U'®V'- V'®U'). (42) 

Thus, Ej = Ei = O for i = 1, 2, . . . , 4 and j = 1, 2, . . . , 8, and the multiplication 
table is^ 

Ml M2 M3 M4 M5 Mg M7 Mg 




The algebra of the Mj's is equivalent to that of 2 x 2 complex matrices [T^ . 

We note that Mi = F and M2 = G, and hence the decomposition for the identity 
tensor is 

I = El + E2 + Ml + M2. (44) 

The Euclidean lengths of the new tensors are || Mi|| = || R2II = V^, || R3 + M4II = || R5 + 
Mg II =2. 

A general symmetric fourth-order tensor of trigonal symmetry is given by the linear 
combination 

A = aEi + 6E2 + c(E3+ E4) + j>Mi + gM2 + r( R3 + M4) + s(M5+ Rg). (45) 
4.7 Rhombic, Monoclinic and Triclinic systems 

The algebra for orthorhombic (equivalently orthotropic) symmetry is the same as that of 
3x3 real matrices plus three real numbers. For monoclinic symmetry it is the same as 
the algebra of pairs of 4 x 4 and 2x2 real matrices. The lowest symmetry, triclinic or no 
symmetry, has algebra the same as that of 6 x 6 real matrices. These low symmetries are 
probably better defined in terms of their group properties, which also serves as efficient 
means for Euclidean projection, as described in Sectional 

^The multiplication table in |19| contains some typographical errors; specifically the elements (3,8), 
(4,5), (5,8), (6,3), (8,4) and (8,6) in the original need to be multiphed by —1. 
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4.8 Six-dimensional representation 

An equivalent 6-dimensional matrix format is described for the various elementary tensors. 
For each symmetry class the 6-dimensional tensors possess the same algebraic properties 
as the fourth-order tensors, with the same multiplication tables and Euclidean lengths. 

The 6-dimensional vectors {a, b, c} represent the second-order tensors {a ® a, b ® b, 
c ® c} that correspond to the axes {a, b, c}, where according to eq. (0), 

a = oiBi + 0262 + 0363 ^ a = (^al, al, al, V2a2a3, V2a3ai, V2aia2j , etc. (46) 
The associated 6x6 base matrices are given next. 

4.8.1 Isotropic system 

The six-dimensional analogs of I, J and IK are I, the 6x6 identity matrix, and J, K, 
where 

J = -(ei + e2 + e3)(ei + e2 + e3)*, K = I-J. (47) 

4.8.2 Cubic system 

The analogs of L and M are 

L = I - aa* - bb* - cc*, M = K - L. (48) 

4.8.3 Transversely isotropic system 

Define the 6— vectors 

P = ^' ^=72*^^ + ^)' ^=72^^"^)' 

w = ^v^ai6i, V2a2^2, aV^^bs, (0263 + 0362), {ash + aibs), (01^2 + 02^ 

then the 6-dimensional matrices for the tensors Ei, E2, E3, E4, F and G are 

El = pp*, E2 = qq*, E3 = pq*, E4 = qp*, 
F = WW* + zz*, G = L — WW*. 

4.8.4 Tetragonal system 

The six-dimensional representation is as for TI, with the addition 

Fi = WW* , F2 = zz* F3 = wz* , F4 = zw* . (51) 

The trigonal system can be treated in the same manner using equations similar to 
for a and (jl^ for w, but the details are omitted for brevity. 



'1, 



(49) 



(50) 
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5 Euclidean projection 

Euclidean projection is an essential ingredient for determining the closest tensors in both 
the Euclidean and the log-Euclidean metrics. It is defined in abstract terms in the next 
subsection, with the remainder of the section devoted to explicit procedures for the pro- 
jection. 



5.1 Definition of the projection 

We wish to find the tensor Csym of a specific symmetry class which minimizes the Eu- 
clidean distance || C — Csymll of an elasticity tensor of arbitrary symmetry, C, from the 
particular symmetry. The solution is a Euclidean decomposition 



-sym 



+ Cxsym , (52) 



where Csym possesses the symmetries appropriate to the symmetry class considered. The 
complement, or residue [3], is orthogonal to Csym, 

( Csym, C_Lsym) = , (53) 

and hence 

II Cf = ||Csym|r+ llCxsymf , (54a) 

II C Csymll II C_|_sym|| • (54b) 

We illustrate the recursive nature of the Euclidean projection for the special case of 
isotropy. Any elasticity may be partitioned into orthogonal components as 

C Cjso ~t~ C_|_iso 

Ciso ~l" Ccub/iso ~l~ Ctct/cub ~l~ Cort/tot ~l~ Cmon/ort ~l~ C_|_mon, (^^) 

where 

CsymB/A = CsymB CgymA, (^6) 

with sym A of higher symmetry than sym B. Equation (j55p partitions C using the path 
from no symmetry to isotropy via cubic symmetry, see Figure ^ An alternative decom- 
position via hexagonal, gives 

C Cjso ~l" Chex/iso ~l~ Ctet/iiex ~l" Cgrt/tct ~l~ Cmon/ort ~l~ C_|_mon' (^'^) 

The different paths in Figure Q were identified by Gazis et al. jS] and discussed more 
recently in greater detail by Chadwick et al. [IHj. The idea of sequential decomposition 
in the Euclidean norm is not new, e.g., t>ut the explicit projection and complement 

have not been represented previously. 

Details of the recursive Euclidean partition of C are given in the Apendix. For in- 
stance, eq. ()A.30|1 represents || C||^ using three different routes from triclinic to isotropic. 
In particular, the length of the isotropic projection is given by 

llCisof = 9/€2 + 20/i2, (58) 

where k and /i are defined in (j23|l . and the distance of C from isotropy is given by 

II Cxisof = II cf - - 20/^2 . (59) 
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It appears that Fedorov |4j was the first to show that the distance from isotropy is mini- 
mized if the isotropic tensor is as given in ()22j). 

The remainder of this Section details three alternative projection methods using basis 
tensors, using symmetry planes or groups, and using 21-dimensional vectors. The pro- 
jection operators are valid for all symmetries, but the method based on symmetry planes 
becomes more complicated for the higher symmetries. Conversely, the projection operator 
using basis tensors becomes simpler at the higher symmetries. The 21-dimensional vector 
procedure is straightforward but is coordinate dependent. 



5.2 Projection using a tensor basis 

We assume that the set of tensors { Vj G Ela} form a linearly independent basis for the 
symmetry sym in the sense that any elasticity tensor of that symmetry may be expressed 
uniquely in terms of linearly independent tensors Vi, V2, . . . Vat, where 2 < < 13 
is the dimension of the space for the material symmetry. We have seen the explicit form 
of several sets of basis tensors in Section |3J e.g., N = 2 for isotropic elasticity. = 13 
corresponds to monoclinic, which is the lowest symmetry apart from triclinic (technically 
= 21) which is no symmetry. The A^ elements of the orthogonal basis are assumed to 
be independent of the elasticity itself, i.e., they do not require "elasticity distributors", 
which are necessary for the spectral decomposition of Ela ^^l- They depend only on 
the choice of the planes and/or axes which define the group G of symmetry preserving 
transformations (see below). The precise form of the basis tensors is irrelevant, all that 
is required is that they be linearly independent in the sense of elements of a vector space, 
and consequently any tensor with the desired symmetry can be expressed 

N 

Csym = J]aiVi. (60) 

i=l 

Minimizing d|,(C, Cgym) with respect to the coefficients a, implies 
_d_ 

Hence, 



^ II C - Cy^f = 2( C, V,) - 2( Csy^, V,) = 0. (61) 



N 

^a,(V„V,) = (C,V,). (62) 

Let D be the N x N symmetric matrix with elements 

A, = (V„V,). (63) 
This is invertible by virtue of the linear independence of the basis tensors, and so 

N 
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Noting that C_Lsym = C — Cgym, we have 

N N 
{ C_Lsym, Csym) = ( C - ^ ttj V,, ^ ttj Yj) 

i=l j=l 
N N N 

= J]a, (C, V,) -Y^Y^a^a, (V„ V,) = 0. (65) 

j=l i=l j=l 

Hence, the partition (jHUjl satisfies the fundamental projection property (jS2I)- It is a simple 
exercise to show that Csym of eq. (j6(J|) minimizes the Euclidean distance || C — Csymll iff 
the coefficients are given by eq. (j64|l . 

The Euclidean projection can also be expressed in terms of tensors of order eight, 
Psym, such that 

Csym = -Psym C, ^ijkl = ^Jklmnrs^mnrs- (66) 

The projector can be expressed in terms of dyadics of basis tensors, 

N 

Psym = Y ® ■ (67) 

We now apply these ideas to particular symmetries, focusing only on the higher sym- 
metries since the lower ones (monoclinic, orthorhombic) are relatively trivial and can be 
handled using the projector of eq. (ffT|l . 

5.3 Projection operators for particular symmetries 

The projector for the symmetries described in Section|4]follow naturally from the vectorial 
formulation. In particular, we note that the D is diagonal in each case, and 

Piso = J ® J + I K ® K, (68a) 
Pcub=JI®J+|L®L+|M®M, (68b) 

Phcx = El ® El + E2 ® E2 + |( Eg + E4) ® ( Eg + E4) 

+ |E(g)E+|G(g)G, (68c) 
Ptet = El ® El + E2 ® E2 + i ( E3 + E4) ® ( E3 + E4) 

+ El ® El + E2 ® E2 + |(E3 + E4) ® (F3 + E4) + i G ® G, (68d) 
Ptrig = El ® El + E2 ® E2 + |( E3 + E4) ® ( E3 + E4) + I Ml ® Ml 

+ iM2® M2 + i(M3+ M4)®(M3+ M4) + i(M5+ M6)®(M5+ Me). (68e) 

Expressions can be given for the lower symmetries, but as discussed, the number of elemen- 
tary tensors involved is much larger. The form of the projectors in ()68|) are independent 
of the crystal axes, and provide a coordinate-free scheme for Euclidean projection. They 
can be easily programmed using the 6-dimensional matrix representation. 

If sym A > sym B, we may convert PsymB into a form that yields CgymB and the 
complement CsymB/A of f|^- A general procedure for achieving this uses the fact that 
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basis tensors for a given symmetry are not unique, and accordingly, alternate forms of 
Pgym may be found. The decomposition 

-PsymB = -PsymA + -PsymB/A; (69) 

can be found by a Gram-Schmid type of process. As an example, let sym A= iso and sym 
B= cub, and consider the tensor basis for cubic symmetry { Vi, V2, V3} = { J, K, L — 
aM}, where a is a free parameter. Requiring that D is diagonal implies a = 3/2, and 
consequently we obtain 

Pcub = J ® J + I K ® K + ^(2L - 3M) ® (2L - 3M) = Piso + ^cub/iso • (70) 

The analogous partition for isotropy and transverse isotropy is described in Appendix A, 
with Phex/iso given in eq. ()A.24jl. 

5.4 Projection using the transformation group 

Let G be the group of transformations for the material symmetry. G C SO (3), the 
space of 3-dimensional special (or proper) orthogonal transformations. Under the change 
of basis associated with Q G SO (3), an elasticity tensor transforms C ^ C' where 
C'ijki = QipQjqQkrQisGpqrs- The projcctiou is 

n 

^sym _ J- \ ^ n^m) ^(m) ^(m) ^{m) ^ x 
ijkl ~ ~ ^ip ^jq ^kr ^Is ^pqrs, I ' -LJ 

m=l 

where n is the order of the group G, and Q^™-*, are elements of the group. The order 
of the groups are as follows (see Table 6 of jSH]): triclinic 1, monoclinic 2, trigonal 6, 
orthorhombic 4, tetragonal 8, transverse isotropy 00 + 2, and isotropy 00^. The projection 
defined by (f7T|) was discussed by, among others, Gazis et al. jHI, and has been applied to 
ultrasonic data by Frangois et al. |T7] . 

The six-dimensional version of the projection is 

n 

Csym = - V Q('")CQ(™)\ (72) 

m=l 

where Q^™) G 5*0(6) is the orthogonal second-order tensor corresponding to Q^™\ satis- 
fying QQ* = Q*Q = I. A necessary and sufficient condition that Cgym has the correct 

symmetry is that Q'^-'-'CsymQ*-"'^ , corresponding to the transformation Q*--'^ is also a mem- 
ber of the symmetry class. Thus, for j = 1, 2, . . . , ra, 

Q(^')C.y^Q(^-)* = - V Q(^')Q(^'CQ« Q(^-)' = - V Q^'^^CQ^'^)' . (73) 
n ^-^ n ^-^ 

i=l k=l 

The latter identity follows from Q^-'^Q*^*^ = Q^'^^ and the fact that the Q*-*^ form a group of 
order n. This proves that Cgym defined in eq. fl71|) is indeed an element of the symmetry 
class. 

The projection formula (f7T|) or (f?^ may be effected in several ways, e.g., by referring 
to the extensive lists of the 6x6 matrices Q*-*-* provided by Cowin and Mehrabadi |35j . 
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We note that each Q*^*^ can be decomposed into a product of reflection operators 
The fundamental operator is R(n) = I — 2n ® n, where the unit vector n is normal to 
the plane. Chadwick et al. [TH] provide a complete description of these sets and groups. 
It is also possible to represent the group of transformations in terms of elements each 
corresponding to a single rotation ITH] . The equivalence between reflection operator 
and rotations follows from the identity Q(e,7r) = — R(e), where Q(n, 6*) G 5*0(3) defines 
rotation about n by angle 9. 

We illustrate this general procedure with an example. 



5.4.1 Example: Projection onto monoclinic symmetry 



There is only one plane of symmetry [18], with normal c, say, and n = 2 group elements 
G = {I, -R(c)}. Using (HH) we have 



mon 1 



Cijki + \Rip{c)Rjg{c)Rkr{c)Ris{c)Cj 



ijkl ~ 2 2 



pqrs ■ 



(74) 



For instance, let c = ea, and using the explicit form of R(e3) from eq. (54) of (where 
it is called R^^^), 

R(e3) = diag(l,l, 1,-1, -1,1), (75) 

then (f?^ gives 

/cii ci2 ci3 cieX 



C^on = ^C + ^R(e3)CR(e3) 



C22 C23 C26 



C33 C36 



S Y M 



C44 C45 



C55 



V 



C66/ 



(76) 



The example illustrates that the group projection method is practical if the order n 
of the group G is not large. This is the case for monoclinic, orthorhombic and perhaps 
trigonal systems, but even for these symmetries it is simpler to implement the algorithm 
by machine. It is not practical for isotropy and transverse isotropy, each with an infinity 
of symmetry planes. Application of (j71|) for cubic isotropy is unwieldy although has been 
performed |17] . The explicit orthogonal basis method of (|6U|) and (|64p provides a simpler 
procedure for the higher symmetries. 



5.5 Projections using 21-dimensional vectors 

The projection is simplest for this representation of Ela because the elastic moduli are 
explicitly represented as a vector. The operator is a 21 x 21 matrix Pgym, and the closest 
elasticity is 

^sym Psyra-^i C^"^) 
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where X is defined in (Q. The various projection matrices can be read off from the 
formulas of Appendix A. Thus, the projector for cubic symmetry is 

|||oooooo 

" i 



cub 



Pcub 09x12 
>0l2x9 0i2xl2, 



Pcub 



? ? ? 



3 3 3 





3 3 3 





yo 



i i i 

* * * 



111 

! ! ! 

! ! ! 

3 3 3/ 



(7J 



Browaeys and Chevrot fl] derived Psym for isotropy, transverse isotropy (hexagonal sym- 
metry), tetragonal, orthorhombic and monoclinic symmetry. They found that in all cases 
except monoclinic that the non-zero part of Pgym reduces to a 9 x 9 matrix Psym, as in 
(jZH))- However , it sh ould be noted that this is not true for the tetragonal projector derived 
here, see eq. (IA.6|1, for reasons discussed in Section 14.51 Nor it is true for the trigonal 
projector, see eq. (jA.lOjl . 



6 Exponential, logarithm and square root of sym- 
metric fourth-order tensors 

Unlike the Euclidean projection for the Frobenius norm, finding minimizers for the loga- 
rithmic and Riemannian distance functions involves evaluation of analytic functions with 
symmetric fourth-order tensor arguments. In this section we derive the necessary ex- 
pressions with emphasis on exponential, logarithm and square root functions that will be 
used to obtain the closest elastic tensors in section [3 For this purpose, let h denote the 
analytic function given by its power series 

oo 

h{x) = amX^- 

m=0 

The higher symmetries are simpler, and are considered first. 



6.1 Isotropic system 

Using (1201) and (ED), for ^ defined in ^ we have 



oo 



h{K) = ^a™(aJ + &K)"^ = ao I + J] a™(a"^ JJ + 6"^ K) 

m=0 m=l 

oo 

= aoI+ ^a^(a™J + 6™K) -ao(JJ+ K) = /i(a) JJ + K. (79) 

m=0 

It follows that 

exp A = J + K. (80) 
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When a 7^ and b ^ 0, the symmetric tensor A is invertible and its inverse is 

A-^ = -J+^K. (81) 
a 

The tensor A is positive definite if a > and b > 0, with square root given by 

A^/2 = v^J + v^K, (82) 

and logarithm 

log A = lnaJ + ln6K. (83) 

6.2 Cubic system 

Using (j2SI)i and for A defined in we have 



oo 



/i( A) = ^ ao(a J + 6 L + c M)™ = ao I + ^ a„(a'" I + b"" iL + c"" M) 



m=0 m=l 

oo 



= aoI+ ^a,r^(a"'J + &™L + c™M) -ao(J+ L+ M) = h{a) J + h{b)l. + h{c)M. 

m=0 

(84) 

In particular, 

exp(A) = e'' J + e^'L + e'^M. (85) 

When a 7^ 0, 6 7^ and c 7^ 0, the symmetric tensor A is invertible and its inverse is 
given by 

A^i = - J+ iL + -M. (86) 

a b c 

If a > 0, 6 > and c > then A is positive definite. Its square root is 

A^/2 = v^J + v^L + VcM, (87) 

and its logarithm is given by 

Log A = lnaJ + ln6L + lncM. (88) 

6.3 Transversely isotropic system 

Using (jn21) and (jH!)), for A defined in (jH^ we have 

00 

h{K) = [aEi + 6E2 + c(E3+ E4)+/F + ^G]™ 

rr!=0 

00 

= ao I + ^ a™ [(a El + 6 E2 + c( E3 + E4))'" + /™' F + g"^ G] 



m=l 

00 



^ a„(a El + 6E2 + c( E3 + E4))™ + h{f) F + h{g) G - ao( F + 



m=0 



/i(0)(Ei+ E2- I) + /i(aEi + 6E2 + c(E3+ E4)) + /i(/) F + /i(^) G. (89) 
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Contrary to the two previous classes of material symmetries, the upper left 4x4 multipli- 
cation table for the elementary tensors of a transversely isotropic medium is not diagonal. 
This fact makes the algebra of transversely isotropic tensors a bit more difficult. Fortu- 
nately, Walpole jinj showed that the algebra for Ei, E2 E3 and E4 is equivalent to the 
algebra of 2 x 2 matrices. The details are described in Appendix El 

For the sake of simplicity of notation, we use the symbol a to denote the quadruple 
(a, b, c, d) and introduce the functions 

(3(a) = -(a-b), 

^ ' 2^ ^' (90) 
7(a) = ^V(a-6)2 + 4(c2 + rf2), 

(5(a) = Va6 - - d? {li ab - - >Q). 

The variable d is not needed in this subsection, so we take d = 0. We used it to keep the 
number of auxiliary functions to a minimum. 

The procedure is to diagonalize the compound tensor argument in 



B = aEi + 6E2 + c(E3+ E4) = (a(a) + 7(a)) E+ + (a(a) - 7(a)) E_, (91) 

where 

E± = -(Ei+ E2) ± -[cosV'(a)(Ei - E2) + sin^^(a)(E3 + E4)], (92) 
with ^/^(a) defined by 

/5(a) c 
cos'?/'(a) = —r^, sin?/^(a) = — — if 7(a) > 0, 
7(a) 7(a) 

^V(a)=0, if 7(a) = 0. 

The derivation is apparent from the results below, and explained in detail in Appendix iBl 
The central feature of this decomposition is that E-t satisfy E-tE-t = E-t, E+ E_ = O, 
E_|_ + E_ = El + E2, and consequently any function, /i, of the tensor B may be expressed 

h{ B) = /i(a(a) + 7(a)) E+ + /i(a(a) - 7(a)) E_ + /i(0)( I - Ei - E2) , (93) 
and hence 

h{ A) = /i(a(a) + 7(a)) E+ + /i(a(a) - 7(a)) E_ + h[j) F + h{g) G . (94) 
Converting back to the Ej's gives 

/i(A) = /i+(a)(Ei+ E2) + /i-(a)[/3(a)(Ei- E2) + c(E3+ E4)] + /i(/) F + /i((7) G , (95) 
where 

^+ ^ /i(a(a)+ 7(a)) + /i(a(a)- 7(a)) 

and 

/z(a(a) -f 7(a)) -/^(a(a)- 7(a)) . ^ ^ . 

h (a) = <; 27(a) 

>'(a(a)), if 7(a) = 0. 
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The results of Appendix El and above imply 

exp(A) = ei(a)Ei + C2(a)E2 + e3(a)(E3+ E4) + e-^ F + G, (96) 

where 

Ci(a) = exp(Q!(a))[cosh(7(a)) + /9(a) sinhc(7(a))], 

e2(a) = exp(a(a))[cosh(7(a)) - /5(a) sinhc(7(a))], (97) 

e3(a) = cexp(a(a)) sinhc(7(a)), 

and sinhc(-) is the hyperbolic sine cardinal function defined by 

fl ifx = 0, 

sinhc(x) = < sinh(x) , (98) 
otherwise. 

Note that sinhc(x) is continuous at 0. 

When a6 — 7^ 0, / 7^ and g ^ 0, the symmetric tensor A is invertible and 

A-^= J &Ei + aE2 -c(E3+ E4)] + ^F + -G. (99) 

ah — J g 

A is positive definite if a > 0, 6 > 0, a6 — > 0, / > and g > 0, with square root 

Ai/2 = \ [{a + (5(a)) Ei + (6 + (5(a)) E2 + c( E3 + E4)] + y^¥ + ^G, 

V2(a(a) + 5(a)) 

(100) 

and logarithm 

LogA= [i(a)Ei + [2(a)E2 + [3(a)(E3+ E4) + In/F + In^G, (101) 

where 

[i(a) =ln5(a)+/3(a)£(a), 

[2(a) = In (5(a) -/5(a)£(a), (102) 
h{a)=ci{a), 

with the function £{■) defined by 

^ if 7(a) = 0, 



tt a 



— -— In ^[^1 otherwise. 

27(a) a(a)-7(a) 

Note that for 7(a) we have 75-^ In ^4^7^-4^ = -rrtanh""^ and that the limit of this 

' ^ ' 27(a) a(a)— 7(a) 7(a) a(a) 

expression is l/a{a) as 7(a) goes to zero. 
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6.4 Tetragonal system 

Using the decomposition (jH^ and the multiphcation table (jHUj) . an analysis similar to 
that used for the transversely isotropic system implies that for a tetragonal tensor A as 
defined in (n?7jl . 

h{A) = /i+(a)(Ei+ E2) +/i-(a) [/3(a)(Ei - E2) + c(E3+ E4)] 

+ /i+(p)(Ei+ E2) + /i"(p)[/?(p)(Ei- E2)+r(E3+ E4)]+%)G. (104) 

As before, a = (a, b, c, d) (with d = 0) and p = (p, g, r, s) (with s = 0). 
Thus, 

exp(A) = Ci(a)Ei + e2(a)E2 + e3(a)(E3+ E4) 

+ ei(p)Ei + e2(p)E2 + e3(p)(E3+ ¥i) + e3Q. (105) 

The tensor A is invertible if afe — 7^ 0, — 7^ and 7^ 0, and 

1 „ _ „ 1 . „ „ „ 1 „ 



-1 



[6Ei + aE2-c(E3+ E4)] + [gFi + J9E2 - r(F3 + F4) 



ao — pq — r'^ g 

(106) 

A is positive definite if a > 0, b > 0, ab — > 0, p > 0, q > 0, pq — r"^ > and (7 > 0, 
with square root 

= /nf f\ ^^ + El + (6 + 5(a)) E2 + c( E3 + E4)] 
V2(a(a) + 5(a)) 

+ , ^ ^1 [(p + 5(p)) El + (g + 5(p)) E2 + r( E3 + E4)] + y^G, (107) 

V2(a(p) + <5(p)) 



and logarithm 



Log A= [i(a)Ei + [2(a)E2 + [3(a)(E3+ E4) 

+ ^l(p) Fi + [2(p) E2 + [3(p)( E3 + F4) + In (7 G. (108) 



6.5 Trigonal system 

The decomposition ()44|) and the multiplication table yield for A defined in fl45p 

h{k) = /i+(a)(Ei + E2) + /i-(a) [/9(a)(Ei - E2) + c(E3+ E4)] 

+ /i+(p)(Ri+ M2) + /i"(p) [/5(p)(Mi - M2) + r(M3+ M4) + s(R5+ Kg)] , 

(109) 

where, as before, a = (a, 6, c, c?) (with d = 0) and p = {p, q, r, s) (here s need not be zero). 
Using these results we have 

exp(A) = ei(a)Ei + e2(a)E2 + e3(a)(E3+ E4) 

+ ei(p)Mi + e2(p)M2 + e3(p)(M3+ M4) + e4(p)(M5+ Me), (HO) 
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where ti{-), 2 = 1,3 are as defined in ^7\i and e4(p) = s exp(a(p)) sinhc(7(p)). The tensor 
A is invertible if ab — ^ and pq — r"^ — ^ 0, with 

A"^ = J J &Ei + aE2-c(E3+ E4)] 
ab — & 

+ \ ^[gMi+pM2-r(M3+ R4)-s(M5+ Me)]. (HI) 

pq — r'' — 

A is positive definite if a > 0, 6 > 0, a6 — > 0, p > 0, g > and pq — r"^ — > 0, with 
square root 

A^/^ = ,^ ^\ ^, = [(a + 5(a)) El + (6 + 5(a)) E^ + c( E3 + E4)] 
v/2(a(a) + 5(a)) 

+ — ==^=== [(p + 5(p)) Ml + (g + (5(p)) M2 + r( M3 + M4) + s( Ms + Mg)] , 
V2(a(p) + 5(p)) 

(112) 

and logarithm 

Log A = [1 (a) El + l2(a)E2 + [3(a)(E3+ E4) 

+ li(p)Mi + l2(p)M2 + [3(p)(M3+ M4) + [4(p)(M5+ Me), (113) 

where [,(■), = 1,3 defined in (fTTl^ and [4(p) = s£(p). 

6.6 Rhombic, monoclinic and triclinic systems 

There are no practical analytical results for the exponential, logarithm and square root 
for these low symmetries. The purely numerical route is recommended. 



7 The closest tensors using logarithmic norms 

We now consider the problem of finding the closest elasticity tensors using the Riemannian 
and the logarithmic distance functions. The solutions use the machinery developed in the 
previous Section for evaluating functions of tensors. We begin with the higher symmetries. 



7.1 The closest isotropic tensor using the Riemannian norm 



We want to find the closest isotropic tensor = 3fi;/j J + 2/i^K to a given elasticity 
tensor C, i.e., find k/j > and fj,R> such that the Riemannian distance dij(C, C^°) is 
minimized, where d/j is defined in eq. ()llcp . To find the optimality conditions, we note 
that minimizing the Riemannian distance is equivalent to minimizing the square of this 
distance, and recall the following result [211 Prop. 2.1] 



d 
It 



tr [Log2 X(t)] = 2 tr 



LogX(t)X-i(t)|x(t) 



Here X(t) is a real matrix- valued function of the real variable t such that, for all t in 
its domain, X(t) is an invertible matrix which does not have eigenvalues on the closed 



MOAKHER AND NORRIS 



February 4, 2008 



25 



negative real line. This result applies to fourth-order elasticity tensors via the isomorphism 
© and ©. 

Using this general result implies two conditions, 

|4( c, Cf) = — tr [Log( C-^ C%') J] , (114a) 



( C, e^°) = — tr [Log( e|°) K] . (114b) 

Hence, the optimality conditions for are 

tr [Log( e|°) J] = 0, tr [Log( C"^ C^°) K] = 0. (115) 
In view of the partition (|2(Jj). addition of these two equations yields 

tr [Log(C-^C^°)] =0, 
which can be rewritten as^ (see Appendix IB]) 

det(c-ie^°) = 1, 

or, equivalent ly, 

det C = det C^°. (116) 

This is the counterpart of the fact that, for the Euclidean distance, the closest isotropic 
tensor has the same trace as the given elasticity tensor. 

7.1.1 The closest isotropic tensor to a given cubic tensor 

For C = aJI + 6L + cM, the conditions (jll5b ) and (jll6|) . after some algebra, reduce to 

In — ^ = 0, ab^c^ = 3/^R(2^R)^ (117) 
a 

which can readily be solved to yield 

3kr = a, 2fin = ip^^f/^ 

7.1.2 The closest isotropic tensor to a given transversely isotropic tensor 

For C = aEi + 6E2 + c(E3 + E4) + / F + ^G, the conditions (fTTHk ) and (ITTH|) yield 

31n5(x)-£(x)(/5(x)-2V2^) = 0, (a6 - c2)/V = 3/^R(2/i^)^ (118) 
where x = (x, y, 2, 0) with 



■^Here and throughout, the determinant of a fourth-order elasticity tensor C is det C = n/=i^J' 
where the A/'s are the eigenvalues of C. From © we have det C = detC, where C is the associated 
six-dimensional second-order tensor. 
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where 

Since S is an orthogonal transformation, it follows that AB — C^ = ab — c^, A + B = a + b, 
xy- z'^ = {ab - c^)/{6krHr) and x + y = {a + 6)/(6/t/j/i/j). 

Equation pi8|) o can be used to eliminate kr in equation piS^ i and hence we obtain a 
single nonlinear equation for a single unknown. An appropriate choice for this unknown 
is the positive variable ^ defined by the ratio 



2fiR 1 - 2ur 

where ur is the Poisson's ratio of the closest isotropic tensor with respect to the Rieman- 
nian distance. Then ()118p yields the following equation for ^ 

^ A + B^^-R{^) AB-C^ ' ^ ^ 

where 



Rio = \j{A-Bey+^c^e, 

and A = A/{det Cy/^, B = B/{det Cy/^, C = C/(det Cy/'^. Note here that for this 
type of symmetry det C = {ab — c^)pg'^. 

The isotropic moduli therefore ^epend upon three unidimensional combinations of the 
transversely isotropic moduli, A, B and C. Once the solution of eq. ^VZ1\^ is found, the 
isotropic moduli are given by 

?>kr = (det C)^/^e^/^ 2fiR = (det Cy/^C^^^ (122) 

7.1.3 The closest isotropic tensor to a given tetragonal tensor 

For C = aEi + 6E2 + c(E3+ E4) +pFi + gFa + r( Fg + F4)+^G, the conditions (fim ,! 
and (fTTEIl yield 

31n(5(x) -£(x)(/3(x) -2V22) = 0, {ab - c^){pq - r^)g^ = 3KR{2fXRy , (123) 

where again x = {x,y,z,0) is defined by ()119|) . The solution is obtained in similar 
manner to that for transverse isotropy, thus, kr and fiR are given by ()122|) but in this 

case det C = {ab — c^){pq — r'^)g'^- 

7.1.4 The closest isotropic tensor to a given trigonal tensor 

For C = aEi + 6E2 + c(E3+ E4)+pMi + gM2 + r(M3+ M4) + s(M5+ Me), the conditions 
(fTTHk ) and yield 

31n5(x) -^(x)(/3(x) -2^2^) = 0, ^{ab - c^){pq - - s^y = 3KR{2fiRf , (124) 

with X = {x,y,z,0) of eq. ()119|) . The isotropic elastic moduli are obtained in similar 
fashion as for of transverse isotropy, i.e., kr and fiR are given by ()122|) but now det C = 
^{ab — c^){pq — — s^)^. 



MOAKHER AND NORRIS 



February 4, 2008 



27 



7.2 The closest cubic tensor using the Riemannian norm 

The optimality conditions that define the closest cubic tensor C^*^ = S/t/? J + 2/iijL + 
2?7/j M to a given elasticity tensor C of lower symmetry, are obtained by minimizing the 
Riemannian distance dij(C, C^^*^): 

C, CS"^) = — tr [Log( Q"") J] = 0, (125a) 

c, CS^b) = —ti [Log( Ct) ^ = 0, (125b) 

^( C, C5j"^) = — tr [Log( C't) M] = 0. (125c) 
oVr Vr 

As K/j > 0, /i/j > and rj^ > 0, these conditions are equivalent to 

tr [Log(C-^C5j"'^) V] =0, V= J, L, M. (126) 

In view of the partition ()26Pi . combining these three equations yields, in the same manner 
as for (P7|l . 

det C = det = 3KR(2/XR)^(2r//j)l (127) 

If C = aEi + 6E2 + c(E3+ E4) +pFi + gFa + r(F3 + F4) +^G, i.e., is of tetragonal 
symmetry, then equation ()12fib ) reduces to 

3 In (5(x) - ^(x) (/5(x) - 2^25) = 0, (128) 

where x = (x, y, z, 0) with 

where S is as defined in p2()|) and 

A = diae 



Again, since S is an orthogonal transformation, it follows that AB — C"^ = ab — c^, 
A + B = a + b, xy — P = {ab — c^) / {Qurtir), and x + y = (a + b) / {Qkrvir). Equation 
(|127|) becomes 

{ab - ^){pq - r'')g^ = 3KRi2fXR)'\2rif. (130) 
Let C and a be the nondimensional positive variables defined by 

then ()128j) and (jl26b ) yield the following coupled two equations for ( and a 



A + Be + Q,{0 ^ C'/V-2(det 



1I/3 



{A - BC) In , : Z - QliO In ^ ' = 0, (131a 



A + BC-Q,{C) ^ ab 



c 



p-ga In— — Q2[a)\n = 0. 131b 

p + qa^-Q2{(r) 
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where 

The three moduh of the closest cubic elasticity tensor are then given by 

3kr = C'/V-^(det C)l/^ 2i^R = C^'/V(det C)^/^ 2r]R = C^^^a^^det C)^/^ 



7.3 The closest transversely isotropic tensor using the Rieman- 
nian norm 

We want to find the closest, in the Riemannian metric, transversely isotropic tensor 



C 



hex 
R 



aEi + 6E2 + c( E3 + E4) + f¥ + gG to a given elasticity tensor C of lower 



symmetry, i.e., find a>0, 6>0/>0, 5'>0 and c with ab — > such that the 
Riemannian distance d/j(C, C^°^) is minimized. 
The optimality conditions are 



(9d2 

da ^« ^ 



dc 



ab — 
2 



ab 



tr [Log(C^^C^'^")(6Ei -CE4)] = 0, 
tr [Log( C"^ C^-^") (a E2 - c E3)] = 0, 



(132a) 
(132b) 



c 



hcx\ 
R ) 



c^) 



tr [Log(C~^C^'='^)(-cEi-cE2 + 6E3 + aE4)] =0, (132c) 
tr [Log( C^*^") F] = 0, (132d) 

(132e) 



ab — 6 
2 

7 



-tr [Log(C-^C^'^")G] = 0. 



Asa>0,6>0, a6 — c^>0,/>0 and (7 > 0, these five conditions are equivalent to 

tr [Log(C-^C^^^) V] = 0, V = 6E1 -CE4, aEs -cEa, El + E2, F, G. (133) 
Combining the conditions for Ei + E2, F andG, and partition p2|l implies the constraint 

det C = det C^^". 

7.4 The closest tensors using the log-Euclidean norm 

The closest elasticity tensors according to the log-Euclidean metric follow by minimizing 



d£(C, Csym)- The stationarity condition implies that (Log C — Log Csym' 
(in the Euclidean sense) to the symmetry class, and hence 

Log Csym = -Psym Log C, 

where the Euclidean projector is defined in ()67p. We may therefore write 

Csym = exp (Psym Log C) . 



is orthogonal 



(134) 



(135) 
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This may be evaluated for the particular symmetries using the explicit expressions for 
Pgym in ()68p. combined with the formulas for the logarithm and exponential of elasticity 
tensors in Section IHl We note in particular the Euclidean property, 

d^ (C, CsymA) = d^ (C, CsymB) + d| (CsymB, CsymA) , symA>symB. (136) 

We note the following identity, which is a consequence of the requirement that (Log C— 
Log C*^^™) is orthogonal to each of the basis tensors, 

tr[ V, Log C] = tr[ V, Log QyJ. (137) 

In the remainder of this Section we apply this to the particular cases of isotropy and cubic 
isotropy, for which we derive explicit formulas for the closest moduli. 

7.4.1 The closest isotropic tensor using the log-Euclidean norm 

There are two basis tensors and therefore the general conditions p37|) become, using (jH^ 
with C|° = 3kl J + 2/iL K, 

ln3/tL = tr[JLog C], ln2/ii = i tr[KLog C]. (138) 

Hence, we obtain explicit formulas 

3Ki = exp(tr[JLog C]), 2/1^ = exp (i tr[KLog C]) . (139) 

Using I = JJ + IK and eq. (plijl . it follows that 

3KL{2fiLf = det C. (140) 

7.4.2 The closest isotropic tensor to a given cubic tensor 

For C = a J + 6L + cM, we use dHHl) and K= L + M to simplify the conditions (fnH|l . 
After some algebra we find ^Tj 

3kl = a, 2fiL = [h'c^f'K 

These moduli coincide with those obtained for the Riemannian norm. This is a conse- 
quence of the fact that when two tensors Ci and C2 commute under multiplication we 
have Log( Ci C2 "'^) = Log Ci — Log C2. 

7.4.3 The closest isotropic tensor to a given transversely isotropic tensor 

For C = aEi + 6E2 + c(E3 + E4) +/E + ^G, we use eqs. flTTTT| . (IT!!H|) i and to 

get 

Iw , 2^ fe-a + 4V2c, a + 6 + 27(a) 

hi?,KL = -Mah - c^) ^ In -H" . (141) 

2 ^ ' 127(a) a + 6- 27(a) ^ ' 

Thus, 

b-n+4v/2c 



3«:z. = , , = , (142) 

+ V(a-6)2 + 4c2 y 
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and the shear modulus then follows from ()140|) as 



nl/5 



(a6-c2)/V 



3kl 



(143) 



Note that the equation 1)14111 for kl can be cast in the form of eq. ()122)l i . where ^ satisfies 
an equation similar to (|12H) . 



A + B-R{1) 



AB-C^ 



0. 



(144) 



7.4.4 The closest isotropic tensor to a given tetragonal tensor 

For C = a Ei + 6 E2 + c( Eg + E4) + p¥i + q ¥ 2 + r ( E3 + ¥4) + g G, we again use eqs. 
([MI), ^Mi and the identity (fX32|l i to derive (fr42|l for kl. The shear modulus is 



2/iz 



(a6 — c^){pq — T^)5'^ 



3kl 



1/5 



(145) 



7.4.5 The closest isotropic tensor to a given trigonal tensor 

For C = aEi + feE2 + c(E3+ E4) +pMi + gM2 + r( M3 + M4) + s(M5 + 
for is again recovered, while /i^ follows from 



1 



3kl 



[ah — c?){pq — r"^ — s 



2\2 



1/5 



), the condition 



(146) 



7.4.6 The closest cubic tensor using the log-Euclidean norm 

Proceeding in the same way as for the isotropic case, in this case with three basis tensors, 
C^"^ = 3KiJI + 2/iiL + 2r]L M, we find that the closest cubic tensor has explicit solution 

3Ki = exp(tr[JLogC]), 2/1^ = exp Q tr[LLog C]^ , 2r/L = exp Q tr[MLog C]^ 

(147) 

The moduli satisfy the same determinant constraint as for the Riemannian norm, in this 

case 

3KL(2/Xi)3(2r/i)2 = det C. (148) 



7.4.7 The closest cubic tensor to a given tetragonal tensor 

For C = aEi+6E2 + c(E3+ E4) + pEi + gFg + r( E3 + F4) + ^G, the bulk modulus kl 
is again given by ()142|1 . We use 



L = Fi + 



(149) 



to obtain 



/6 2/3 [ P + q+Vip-qy + ^r ^^ ^vT;^^)^ 



p + q — ^/ {p — qY + 4r2 



(150) 
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The shear modulus rji follows from 



2VL 



{ah — c^){pq — r'^)g'^ 
3«:L(2/ii)3 



1/2 



(151) 



7.4.8 The closest transversely isotropic tensor using the log-Euclidean norm 

In this case there are five basis tensors, C^'^^ = aEi + 6E2 + c(E3 + E4) + f¥ + gQ. 
Using eqs. (jlUlj) and (jl37p . we have 

(152) 



/ = exp ( -tr[FLog 



^ = exp ( -tr[GLog C] 



and the remaining three moduli follow from the identities 

[i(a) =tr[EiLogC], ^(a) = tr[E2 Log C], 13(a) = ^ tr[( E3+ E4) Log C]. (153) 
Thus, 

^cosh(6' + </)cosh6') , ^cosh(6' — <7!)cosh6') ^sinh({/)cosh6') , ,^ 
a = d — , = — , c = d —„ , if 7^ 0, 



cosh 9 

a = Se^, b = Se-^, c = 0, 
where 

= ^tr[(E3+ E4)Log C], 



cosh 6' 



cosh 9 



(154) 



if 6 = 0, 



1 



^ = _tr[(Ei- E2)Log C], 



(5 = exp ( - tr[( El + E2) Log C] 



^ = log 



1 



1/2 



(155a) 
(155b) 



8 Application and numerical examples 

We illustrate the methods developed above by considering an example of a general elas- 
ticity tensor, with no assumed symmetry. This type of data raises a problem typically 
encountered, i.e., find the optimal orientation of the symmetry axes in addition to finding 
the closest elastic tensors for a given symmetry axes or planes. We first present the data, 
and the isotropic approximations, and then consider the question of orientation. 



8.1 Example 

A complete set of 21 elastic constants were determined ultrasonically by Frangois et al. 
[7j. The reader is referred to their paper for details of the measurement technique. The 
raw moduli are 

/243 136 135 22 52 -17\ 
-28 11 16 
29 -49 3 
-28 29 133 -10 -4 
-49 -10 119 -2 
3 -4 -2 130/ 



136 239 137 

135 137 233 
22 

52 11 

\-17 16 



(GPa). 



(156) 
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We first search for tlie presence of symmetry planes, which are an indicator of underlying 
symmetry. Following Cowin and Mehrabadi [33], define A and B by 



If A and B have no common eigenvectors then there are no planes of reflection symmetry 
and the material has no effective elastic symmetry jlHj. We find that the two sets of eigen- 
vectors are not coincident, and the smallest angle between any pair from the two sets of 
eigenvectors is 16°. Materials with symmetry higher than orthorhombic have five or fewer 
distinct eigenvalues A/ [20] • For the given moduli, we find A/ = 47, 79, 244, 285, 312, 512, 
consistent with the material having no symmetry plane. 

We next consider isotropic approximations to the 21 moduli. The Euclidean and log- 
Euclidean approximations follow from eqs. ()68a|) and p39|) . respectively. The numerical 
procedure for finding the closest, in the Riemannian norm, isotropic tensor C^^° = 3k/j J + 
2fifi K is as follows. Define the function 

F{3kr, 2fiu) = tr [Log( C"^ (St^R I + 2fiR K) ) J] . 

By the equal determinant rule pi6|) . we have ?)Kr{2^r)^ = det C, and therefore we need 
only solve the equation for one of kr or fiR. In practice, we solve 

F{3kr, (tf )'^') = 0' 

using Newton's method. 

The isotropic approximations are 









Euclidean, C 


170.11 


96.87 


Euclidean, S 


169.33 


55.81 


log — Euclidean 


169.84 


75.91 


Riemannian 


169.69 


75.92 



'158) 



The Euclidean projection is obviously not invariant under inversion. Another way to see 
this is to consider the product 

SisoCiso ^ 1.001 + 0.74K, (159) 

where Ciso and Siso are the isotropic projections for the stiffness C and compliance S = 
C~^, respectively. The product (jl59|) is an isotropic tensor, as expected, but not the 

identity. The analogous product SsymCsym of the Euclidean stiffness and compliance 
projections becomes closer to the identity as the symmetry is reduced, although not 
uniformly. As one measure, we note that the parameter det (SsymCsym), which is unity for 
the full triclinic matrices, takes the values 4.2, 4.6, 6.6, 4.5, 8.1, 5.3, 15.8 for mon, ort, 
trig, tet, hex, cub and iso, respectively. 



8.2 Orientation effects 



This set of moduli present the general problem of finding the optimal orientation of the 
symmetry axes, for given symmetries. This can be done by a "brute force" approach of 
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searching over all possible orientations of the basis vectors. In practice this is achieved 
using Euler angles (6*1, 6*2, ^3) to transform from {ei, e2, e^} {e'j^, e'g, e'^} by first rotating 
about the 63 axis by ^1, then about the intermediate e'/ axis by ^2, and finally about the 
63 axis by 9^. The moduli transform as C — > C, and for each triple (^1, ^3) the 
closest elasticity tensors Cgym of different symmetries are found. Numerically, we fix the 
symmetries by reference to the original basis, with {a, b, c} = {ei, e2, 63}. 
Let d be any of the three distance functions, and define 

d (C, Ciso) 

Thus, < psym < 1, with equality only if the rotated C' has symmetry sym. The 
preliminary analysis above indicates the absence of any symmetry so we expect Psym to 
be less than unity. We also define 

Plym = max PsyraiOi, 02, 63). (161) 

91,92,03 

Orientations at which psym = Psym are candidates for symmetry axes that best approximate 
the moduli. 



dp, C 


Pcub 


Phcx 


Ptet 


Port 


Ptrig 


Pmon 


cub 


0.91 


0.60 


0.91 


0.92 


0.61 


0.94 


hex 


0.85 


0.60 


0.91 


0.87 


0.62 


0.93 


tet 


0.01 


0.45 


0.92 


0.88 


0.47 


0.94 


ort 


0.03 


0.54 


0.90 


0.94 


0.58 


0.95 


trig 


0.12 


0.20 


0.21 


0.20 


0.95 


0.22 


mon 


0.02 


0.07 


0.35 


0.81 


0.08 


0.98 



Table 1 • Psym calculated using the Euclidean projection of the stiffness C. The diago- 
nal elements in bold are the values of /Ogym for the symmetries indicated in the 
left column. The other numbers in each row are the values of Psym evaluated 
at the same orientation as /sLm- 



Tables 1 and 2 list the results using the Euclidean projection for the stiffness and compli- 
ance, respectively. Table 3 gives the results for the log-Euclidean distance. The results in 
Tables 1-3 were obtained using 60 x 60 x 60 discretized Euler angles. In general, larger 
values give an indication of the proximity of the symmetry. These numbers taken together 
suggests that the material is not well approximated by hex, but that cub and certainly tet 
are reasonable candidates for approximating symmetries. The analogous log-Euclidean 
computations, given in Table 3, reinforce this view. Also, p^^^ = 0.95 for the Riemannian 
norm, using the method below. We focus on the cubic approximation for the remainder 
of this Section. 
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Pcub 


Phcx 


Ptct 


Port 


Ptrig 


Pmon 


cub 


0.82 


0.54 


0.83 


0.97 


0.54 


0.99 


hex 


0.18 


0.71 


0.90 


0.83 


0.72 


0.97 


tet 


0.82 


0.35 


0.96 


0.97 


0.35 


0.98 


ort 


0.82 


0.35 


0.96 


0.97 


0.35 


0.98 


trig 


0.21 


0.32 


0.33 


0.42 


0.84 


0.43 


mon 


0.82 


0.54 


0.83 


0.97 


0.54 


0.99 



Table 2: Psym calculated using the Euclidean projection of the compliance C . The 
elements are determined in the same manner as in Table 1. 



dL 


Pcub 


Phcx 


Ptet 


Port 


Ptrig 


Pmon 


cub 


0.92 


0.68 


0.94 


0.96 


0.69 


0.96 


hex 


0.87 


0.69 


0.93 


0.92 


0.70 


0.96 


tet 


0.11 


0.47 


0.95 


0.71 


0.48 


0.96 


ort 


0.09 


0.20 


0.37 


0.96 


0.21 


0.98 


trig 


0.26 


0.43 


0.43 


0.44 


0.94 


0.47 


mon 


0.91 


0.60 


0.92 


0.96 


0.62 


0.99 



Table 3: Psym calculated using the log-Euclidean distance function, determined in the 
same manner as in Table 1. In particular the diagonal elements are /Ogym- 

8.3 Cubic approximations 

We note some properties of the optimally oriented cubic approximations. First, the de- 
composition ()70|) for the Euclidean cubic projection has an interesting imphcation. Noting 
that 

2L-3M = 5JI + 2K-5(a(g)a(g)a(g)a + b(g)b(g)b(g)b + c(g)C(g)C(g)c), (162) 
and using (|58|) . we may write the cubic length 

II Ceubf = 9/^2 + 20p2 + |(3k + 4/i - Caa - Cfcb - c,,f , (163) 

where k and \i are Fedorov's isotropic moduli, and Caa = { C, a®a®a®a) = Cijkiciiajakai, 
etc. Consider || Ccubll as a function of the orientation of the cube axes. Since k and fi are 
isotropic invariants, and Caa > 0, Cbb > 0, Ccc > on account of the positive definite nature 
of C, it follows that the largest length occurs when Caa + Cbb + Ccc achieves it smallest 
value. This implies that the best cubic approximation in the Euclidean sense occurs in the 
coordinate system with smallest value of {c'n + + C33). 

The closest cubic material for log-Euclidean distance function is that which minimizes 
dL(C, Ccub)- Using the Euclidean property, see ()136|1 . that 

dl ( C, Ceub) = di ( C, Ci,o) - dl ( Ccub, Ciso) , (164) 

and the fact that dii^C, Cjso) is unchanged under rotation, it follows that the optimal 
Ccub maximizes di(Ccub7 Ciso)- Let p^j denote the isotropic modulus, from eq. (jl39p 9. 
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then since the bulk modulus is the same for isotropy and cubic symmetry, it follows 
that 

dL(Ceub, Ci,o) = A/f llog^l (165) 

Therefore, the optimal orientation is that for which + ^) achieves its largest value. 

Alternatively, if we define the Euclidean cubic approximation for C as = 3^^? J + 
+ 2?7^L, then up = k, and the optimal fip and rip maximize Ifip — Vf\ subject to 
the additive constraint + ^Vf = where /i is the isotropic (Fedorov) shear modulus. 
By comparison, the optimal log- Euclidean moduli maximize \ \ogfiL — log?]^! subject to 
the constraint | logfip + | log r]L = log /i. 

The closest cubic tensor in the Riemannian norm, 
is obtained using a numerical scheme similar to that for the isotropic case. Define the 
functions 

FiiS^n, ^fiR, 2r]R) = tr[Log( £-^3^^ J + 2^^L + 27]rM)) J], 
F2{3kr, 2fiR, 2r]R) = tr[Log( C'\3kr^ + 2fiRh + 27]rM)) L]. 

Using the condition p27j) of equality of determinants we can eliminate one of the three 
unknowns, and therefore only need to solve two equations for kr and fiR 

Fii^-R^ 2;.^, (3;^)^/^) = 0, F,{3kr, 2/.^, (3,^^)^/^) = 0. 

These equations are solved by the two-dimensional Newton's method. 

The numerical search for the optimal cubic approximation indicates that the orienta- 
tion which yields the maximum p*^y^ coincides for the three distance functions for this set 
of moduli. The optimal cubic moduli are as follows 







/i 


V 


Eucl, C 


170.1 


139.7 


32.6 


Eucl, § 


169.3 


135.1 


29.7 


log — Eucl 


169.7 


137.5 


31.2 


Riemannian 


169.8 


138.1 


30.9 



(166) 



The Euclidean moduli of the first row are in agreement with Frangois et al. j7j, who used a 
method based on transformation groups to effect the projection. Note that the Euclidean 
projection for compliance again gives a different answer than for stiffness. Also, the bulk 
modulus for the log-Euclidean distance is the same as for the isotropic approximation, see 

(UnHD. 

This example illustrates the practical application of all three distance functions. The 
procedure is similar for materials in which the closest or best symmetry for approximation 
might be, for instance, transversely isotropic. The Euclidean projection follow from eq. 
fl68c|) and the log- Euclidean moduli from eqs. ()152j) and p53p . The closest transversely 
isotropic tensor in the Riemannian norm can be found in a similar manner as above - in 
this case a system of four nonlinear equations for four unknowns is obtained, which can 
be solved by Newton's method. 
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9 Conclusions 

Three distance functions for elastic tensors have been introduced and their apphcation to 
elasticity examined. The Euclidean distance function is the most commonly used, but it 
lacks the property of invariance under inversion. Two distance functions with this impor- 
tant property have been described in detail: the log- Euclidean and the Riemannian norms. 
For each of the three distance functions we have developed a coordinate-independent pro- 
cedure to find the closest elasticity tensor of a given symmetry to a set of moduli of a 
lower symmetry. 

For the Euclidean distance function we have described a projection scheme using basis 
tensors, similar to vector space projections. Explicit forms for the projection operator have 
been given for isotropic, cubic, transversely isotropic, tetragonal and trigonal symmetries. 
The projection method was compared with the group transformation scheme and with 
the 21-dimensional approach. We also described the form of the projection operator using 
6x6 matrices, which is easily implemented on a computer, and derived explicit expressions 
for the tensor complements/residues, and for the lengths of the projections. This allows 
one to decompose, in a Pythagorean sense, the elastic stiffness or compliance, although 
with different results for each. 

Detailed and practical results have been presented for applying the logarithmic based 
distance functions to the same problem of approximating using a higher elastic symmetry. 
Both the log-Euclidean and the Riemannian distance functions are invariant under inver- 
sion. As such they offer an unambiguous method for determining the closest moduli of a 
given symmetry, providing unique results regardless of whether stiffness or compliance is 
considered. 

Calculation of the logarithmic distance functions requires evaluation of the logarithm, 
square root and exponential of tensors. Semi-analytical procedures for obtaining these 
for the most important, higher, elastic symmetries have been derived. For the first time, 
practical expressions are available for computing the logarithmic distances for symme- 
tries higher than orthotropic, that is: isotropic, cubic, hexagonal (transversely isotropic), 
tetragonal, and trigonal. The numerical example considered illustrates how a data set 
with no symmetry can be reduced to find the optimally fitting cubic moduli. The Eu- 
clidean/Frobenius norm yields different results based on whether stiffness or compliances 
is use, but the log-Euclidean and the Riemannian distance functions provide unique an- 
swers. 
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Appendix 

A The closest tensors using the Euclidean norm 

Results for the Euclidean projection are summarized. We first give the explicit form of 
the projection for the various symmetries. 



A.l Summary of the Euclidean projections 

The projections onto the symmetry classes are expressed in terms of the elements of 
the six-dimensional second-order tensor C of eq. ((Tj). For the purpose of subsequent 
calculations, we take {a, b, c} = {61,62,63}. The symmetry classes are characterized 
by a distinct direction, the direction of the monoclinic symmetry plane, the direction 
perpendicular to the plane spanned by the normals to the planes of reflection symmetry 
for trigonal and tetragonal symmetry, the axis of transversely isotropic symmetry, a cube 
axis. In each case we take this direction as 63. 



A. 1.1 Monoclinic and orthorhombic symmetry 
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/cii 


C12 


Cl3 








cw\ 




C22 


C23 








C26 






C33 








C36 
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C45 
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M 




C55 





V 










C66/ 


/cii 


C12 


Cl3 








o\ 




C22 


C23 















C33 

















C44 








S 


Y 


M 




C55 






C66/ 



(A.i: 



(A.2) 



A. 1.2 Tetragonal symmetry 

Ctet as defined in (pTTjl becomes in 6 x 6 notation 



/Ub + q) \{h-q) 



Y 



V2 



c 



V2' 

a 
M 







9 








9 



\ 







p I 



(A.3) 



and the norm is 



+ q^ + 



2/. 



(A.4) 
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Performing the inner products using with N = 7 and { Vi, V2, V3, V4, V5, Ve, Ve} 
{El, E2, (E3 + E4), Fi, F2, (F3 + F4), G} gives 



a = C33, 6 = -(cii + C22 + 2C12) 



P = Cq6, g = -(cii + C22 - 2C12) 



1 

1 



— (C13 + C23), g = ^{ca4 + C55), 



r = -^(Ci6 -C26). 



In summary, the projection is 

/|(C11 + C22) C12 |(C13 + C23) 

|(C11 + C22) |(C13 + C23) 



C 



tct 



C33 







r(c44 + C55) 



(A.5) 











SYM 



v 



^(C44 + C55) 



|(C16 - C26)\ 
|(C26 - Cie) 





C66 

(A.6) 



A. 1.3 Trigonal symmetry 

Ctrig as defined in becomes 

/Ub + P) Ub-P) 0\ 



•'trig 



a 



V2' 





SYM 





(I 





— s 
r 

P / 



(A.7) 



and the norm is 



where 



-tct I 



+ 6^ + 2c^ + + 2g2 + + As' 



1 1 

a = C33, 6 = -(cii + C22 + 2612), c = --^(ci3 + C23), P = cee, 

In this case, 



(Ci5 — C25) — -C46. 



(A.8) 



(A.9) 



Ctrig 



fell 


cIi-cIq |(C13 + C23) 






o\ 




Cn |(C13 + C23) 











C33 















^(C44 + C55) 





^46 




SYM 




i(c44 + C55) 


^56 


\ 








C66/ 



(A.IO) 
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where 



'11 



-66 



-46 



-56 



I (3Cii + 3C22 + 2ci2 + 2C66) , 
I (Cll + C22 - 2ci2 + 2C66) , 



9 I - -^Cis 



9 I C56 + ;^Ci4 



75^25 
75^24 



(A.lla) 
(A.llb) 

(A.llc) 
(A.lld) 



A. 1.4 Transverse isotropy 

With Chox defined as in (jH^ . we have 



f^b + f) \ip-f) 75C 



^hox 



\{h + f) 

a 



9 



SYM 



and 







9 
// 



c 



hex 



6^ + + 2/2 + 



/ = + g) = i(cii + C22 + 2c66 - 2C12) 



Thus, 



(A.12) 



The projection onto the basis tensors { Vi, V2, V3, V4, V5} = { Ei, E2, ( E3 + E4) 
yields a, b, c and g as given in ()A.5jl and 



C 



hex 



-11 "-11 



-66 



-11 



^(Cl3 + C23) 
2 (Cl3 + C23) 
C33 







;(C44 + C55) 

















SYM 



\ 



i(c44 + C55) 

^66/ 



(A.13) 

F, G}, 

(A.14) 



(A.15) 



A. 1.5 Cubic symmetry 

The 6-D matrix associated with Ccub defined in (|28p is 



/Ua + 2c) 



C 



cub 



Ua + 2c) |(a-c) 



SYM 



with length 



■'Cub 



\ 



(a + 2c) 

6 

b 



36^ + 2c2 



(A.16) 



(A.17) 
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The Euclidean projection onto cubic symmetry follows from ()64j) with = 3 and { Vi, V2, V3} 
= { J, L, M}, and gives 



a = |(cii + C22 + C33 + 2ci2 + 2ci3 + 2C23) 
b = I (£44 + C55 + Cge) , 

C = |(Cll + C22 + C33 — C12 — C23 — C31). 

Note that a = 3k where k defined in eq. ()2Hj) i . In summary, the projection is 



■^cub 



where 



'11 



o(Cll + C22 + C33J 



-12 



/5fi 


AC 
•-12 
AC 

•^11 


AC 
•-12 
AC 
•-12 

AC 

•-11 












o\ 







SYM 




AC 

•^66 




AC 

He 
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^66/ 


= 1 


[Cl2 + 


Cl3 + 


C23), 


aC 
'-66 





C44 + C55 + C66j 



(A. 18a) 
(A.18b) 
(A.lSc) 



(A.19) 



(A.20) 



A. 1.6 Isotropy 

With Ciso as given by we have 



SYM 



v 









o\ 




















2/i 



















2fi) 



(A.21: 



where k and /i are defined in These follow from the general projection formula (jMj) 
with iV = 2 and { Vi, V2} = {J, K}. 



A. 2 The space between isotropy and transverse isotropy 

We use a Gram-Schmid approach to define an orthogonal subspace between isotropy and 
transverse isotropy (hexagonal symmetry). Noting that 



\ [Ei + 2E2 + v^(E3+E4)], K=\ [2Ei+E2-y2(E3+E4) + 3F+3( 



we introduce 



Li = G - F, 
L2 = 2 El - 2 E2 



V2 



E3 + E4) 



L3 = 8E1 + 4E2 -4V2(E3 



E4 



3F- 31 



(A.22) 

(A.23a) 
(A.23b) 
(A.23c) 
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Thus, (Li, hj) = 0, i ^ j, and (L^, J) = 0, (L^, K) = 0, i = 1,2,3. Hence {Yi,i = 
1,2,... 5} = {J, K, Li, L2, L3} form an basis for Hex, with D = diag(l, 5, 4, 9, 180). 
Equation ^7\i gives 



Phey, = J®J+^K®K+^Li(8)Li + iL2(g)L2 + -^I^3®L3 
5 4 9 180 

-Piso ~l~ -Phcx/iso • 



(A.24) 



Thus, we have the exphcit decomposition into orthogonal complements Chcx = Cjso 
where 

hex/iso = dihi + d2 IL2 + (is L3, 



Chex/iso, where 



Chex/iso|r = K + 9^2 + 1804 



(A.25) 



with 



di = \{g-f), d2 = l{2a-2b+V2c), d^ = ^{4a + 2b - 4V2c - 3f - 3g), (A.26) 
and a, . . . ,g are given by ()A.5|) and ()A.14|) . Eliminating the latter gives 

rfi = I [2(c44 + C55 + C12) - (cii + £22 + 2c66)] , (A.27a) 
d2 = l [ci3 + C23 + 2c33 - (cn + C22 + 2ci2)] , (A.27b) 
^3 = 3B0 [^11 + ^22 + 14ci2 + 16(c33 - Ci3 - C23) - 6(c44 + C55 + Cee)] . (A. 27c) 

A. 3 Complements 

Some orthogonal complements or residues [3] between the symmetry classes are presented 
here. These lead to the explicit expressions for elastic lengths in the next subsection, eqs. 
flO0|l and (HSH). Thus, 

/ Ci4 Ci5 \ 
C24 C25 
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S Y M C56 

V 0/ 
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S Y M 
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(A.28c) 
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C 
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(A.28d) 
(A.28e) 



\ 
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(A.28f) 



(A.28g) 



where 



ro = |(cii + C22 - 2ci2 - 2c, 
r2 = |(ci3 + C23 - 2C12), 
r = (cii + C22 + C33 - C12 



66 J 



Cl3 



ri = 
^"3 = 
C23 



|(cii + C22 - 2C33), 

|(C44 + C55 - 2Cqq), 

■ C44 — C55 — Cee) • 



(A.29a) 
(A.29b) 
(A.29c) 



Also, Chcx/iso follows from ()A.27|) and ()A.25|) but is not given explicitly because of its 
length, and neither are the complements involving Ctrig- 



A. 4 The Euclidean distances 

Three projections of the elastic length are as follows 
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These can be explicitly calculated using 



^mon/ort 



c 



ort/hcx 
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mon /tot 
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tet / cub 



c 



cub/iso 
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_r) / ^2 I ^2 I ^2 I ^2 I ^2 I ^2 I ^2 i ^2 "N 
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'2 (Ci6 + + + C45) , 



(A.31a) 
(A.31b) 



4 (Cii + C22 — 2Ci2 — 2C66)^ + ^(Cii — £22)^ + (Ci3 — £23)^ + ^(£44 — C55) 



(A.31c) 



4(Cll - £22)^ + (Cl3 - £23)^ + |(C44 - £55)^ + (Ci6 + C26)^ + 2c^5 + 2c: 



36 ' 



and eqs. 

II C- 



|(C11 + C22 - 2C33)2 + |(C13 + C23 - 2ci2)^ 
+ |(C44 + C55 — 2CqqY + (Ci6 — C2gY, 

15 (Cii + C22 + C33 — C12 — Ci3 — C23 — C44 — C55 — Cqq) , 

I (Cll + C22 — 2Ci2 — 2C66)^ + (Ci6 — £25)^ 

■Jq [2(£44 + £55 + £12) — (£11 + £22 + 2£66)]^ 
+ I [£i3 + £23 + 2£33 - (£11 + £22 + 2£i2)]^ 

+ [Cll + £22 + 14£l2 + 16(£33 - £13 - £23) - 6(£44 + £55 + £66) 



and give 

=1 (£11 + £22 + £33 + 2£i2 + 2£i3 + 2£23)^ 

+ ^ [2(£ll + £22 + £33 ~ £12 ~ £23 ~ £31) + 3(£44 + £55 + £66)]^ • 



(A.31d) 

(A.31e) 
(A.31f) 

(A.31g) 



(A.31h) 



(A.32) 



All versions of the elastic length in ()A.30|) are the sum of 21 positive numbers, the same 
as the raw form of the squared length in any rectangular coordinate system using the 
Voigt elements. However, the three sets of 21 positive numbers in flA.30|) contain far more 
information about the underlying elastic symmetry of the material. 



B Exponential, logarithm and square root of Hermi- 
tian matrices 



Let A be the 2x2 Hermitian matrix 

A = 



a c + id 
c — id b 



(B.i: 



where a, b, c and d are real numbers. The matrix can be decomposed into the sum of an 
isotropic part and a deviatoric part. 



A = A, 



iso T ^dcv 



A, 



a 
a 
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(3 c + id 
c — id —(3 



(B.2) 
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where a = (a + b)/2, jS = {a — b)/2 and 7 = \/ (3'^ + + d?. Note that A^cv vanishes if 
and only if 7 = 0, i.e., a = b and c = c? = 0, and that A\^^ = 7^/. The eigenvalues of A 
are real and given by 

A± = a ± 7. 



For convenience we introduce the matrix 

A'-- 



I if 7 = 0, 

^Adev if 7 > 0, 



(B.3) 



[B.4) 



so that A'"^ = I. The exponential of A is then 

( e"/ 

exp A = e"(cosh7/ + sinh A') = < sinh7 

e" (cosh 7/ H 

7 



if 7 = 0, 
Adev) if 7 > 0. 



(B.5) 



When ab — ^ the matrix A is invertible and its inverse is 



A~' 



ab — — d? 
1 

ab — — d? 



b —{c + id) 
-(c — id) a 



ab — — 



[al - 7A') 



[al - A 



dcv J ■ 



(B.6) 



If a > 0, 6 > and ab — — d^ > 0, then A is positive definite. Let 5 = \J ab — (? — rf^, 
then the positive-definite square root of A is 



^1/2 



a + 5 id 
c — id 6 + 5 



and the inverse of its square root is given by 



[(a + 5) / + A 



dcv I 



(B.7) 



^-1/2 ^ ^ 
The logarithm is 



6 + 5 -{c^id) 
-(c — a + (5 



In 5/ 



5v/2(a + <5) 



[(« + 5)J-Adev]. (B.8) 



if 7 = 0, 

ln5/ + -ln J^^^Adev if7>0. 



LogA = In 5/ + In. /^^^ A' 
V q; — 7 

7 V 7 

More generally, for any analytic function / we have 

r/(a)/ if 7 = 0, 

j{A) = l /(c + 7) + /(a-7) ^ /(c, + ^)-/(c-7) ^^ . 

I 2 27 ^ ■ 

The Cayley-Hamilton theorem states that 

det(A)J-tr(A)A + A^ = 0. 



(B.9) 



(B.IO) 



(B.ii: 
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Therefore, any analytic function of A can be written as a linear combination of /, A and 
A"^, or, if A is invertible, as a linear combination of /, A and A~^. Using the fact that 
2A^ev = A — 6'^A~^, we have the alternative expressions for the exponential, logarithm 
and square roots 



exp A 



Log A = In 61 + 

A'/' ^ ^ 

^-1/2 



In 



27 

a + 7 



v/2(« + 5) 
1 



47 a — 7 
[61 + A] , 



{A - 6' A-') 



[I + 6 A-'] . 



V2(a + 6) 

We recall that for any n x n matrix M we have 

det(expM) =e^'^. 
Therefore, when M = Log P for some matrix P we get 



(B.12a) 
(B.12b) 
(B.12c) 

(B.12d) 
(B.13) 



detP = e''^^°^^\ 



(B.14) 
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